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We  theoretically  investigate  the  flow  of  electromagnetic  waves  in  complex  honeycomb  photonic  lattices  with 
local  VT  symmetries.  Such  VT  structure  is  introduced  via  a  judicious  arrangement  of  gain  and  loss  across  the 
honeycomb  lattice,  characterized  by  a  gain  and  loss  parameter  y .  We  found  a  class  of  conical  diffraction  phenom¬ 
ena  where  the  formed  cone  is  brighter  and  travels  along  the  lattice  with  a  transverse  speed  proportional  to  *Jy . 
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I.  INTRODUCTION 

Conical  refraction  phenomena,  i.e.,  the  spreading  into  a 
hollow  cone  of  an  unpolarized  light  beam  entering  a  biaxial 
crystal  along  its  optic  axis,  are  fundamental  in  classical  optics 
and  in  mathematical  physics  [1-4].  Originally  predicted  by 
Hamilton  in  1837  [3]  and  experimentally  observed  by  Lloyd 
[4],  these  phenomena  have  been  intensively  studied  in  recent 
years  by  a  large  community  of  theorists  and  experimentalists 
[1-8].  The  physical  origin  of  the  phenomenon  is  associated 
with  the  existence  of  the  legendary  diabolical  points ,  which 
emerge  along  the  axis  of  intersection  of  the  two  shells 
associated  with  the  wave  surface.  Around  a  diabolical  point 
the  energy  dispersion  relation  is  linear  while  the  direction  of 
the  group  velocity  is  not  uniquely  defined.  Recently,  conical 
diffraction  was  observed  in  two-dimensional  photonic  honey¬ 
comb  lattices  [5]  which  share  key  common  features,  including 
the  existence  of  diabolical  points,  with  the  band  structure  of 
graphene  in  condensed  matter  physics  literature.  In  graphene, 
the  electrons  around  the  diabolical  points  of  the  band  structure 
behave  as  massless  relativistic  fermions ,  thus  resulting  in 
extremely  high  electron  mobility.  Both  photonic  and  electronic 
graphene  structures  allow  us  to  test  experimentally  various 
legendary  predictions  of  relativistic  quantum  mechanics  such 
as  the  Klein  paradox  [9]  and  the  dynamics  of  optical 
tachyons  [10]. 

While  diabolical  points  are  spectral  singularities  associated 
with  Hermitian  systems,  for  pseudo-Hermitian  Hamiltonians, 
like  those  used  for  the  theoretical  description  of  non-Hermitian 
optics,  a  topologically  different  singularity  may  appear:  an 
exceptional  point  (EP),  where  not  only  the  eigenvalues  but 
also  the  associated  eigenstates  coalesce.  Pseudo-Hermitian 
optics  is  a  rapidly  developing  field  which  aims,  via  a  judicious 
design  that  involves  the  combination  of  delicately  balanced 
amplification  and  absorption  regions  together  with  the  mod¬ 
ulation  of  the  index  of  refraction,  to  achieve  new  classes  of 
synthetic  metamaterials  that  can  give  rise  to  altogether  new 
physical  behavior  and  novel  functionality  [11,12].  The  idea 
can  be  carried  out  via  index-guided  geometries  with  special 
antilinear  symmetries.  Adopting  a  Schrodinger  language  that  is 
applicable  in  the  paraxial  approximation,  the  effective  Hamil¬ 
tonian  that  governs  the  optical  beam  evolution  is  non- 
Hermitian  and  commutes  with  the  combined  parity  ( V ) 
and  time  (T)  operator  [13-15].  In  optics,  VT  symmetry 


demands  that  the  complex  refractive  index  obeys  the  condition 
n(r)  =  n*(— r).  It  can  be  shown  that  for  such  structures,  a 
real  propagation  constant  (eigenenergies  in  the  Hamiltonian 
language)  exists  for  some  range  (the  so-called  exact  phase) 
of  the  gain  and  loss  coefficient.  For  larger  values  of  this 
coefficient,  the  system  undergoes  a  spontaneous  symmetry 
breaking ,  corresponding  to  a  transition  from  real  to  complex 
spectra  (the  so-called  broken  phase).  The  phase  transition 
point  shows  all  the  characteristics  of  an  exceptional-point 
singularity. 

VT  symmetries  are  not  only  novel  mathematical  curiosi¬ 
ties.  In  a  series  of  recent  experimental  papers,  VT  dynamics 
were  investigated  and  key  predictions  confirmed  and  demon¬ 
strated  [16-19].  Symmetry  breaking  has  been  experimentally 
observed  in  non-Hermitian  structures  [16-18],  while  power- 
law  growth — characteristic  of  phase  transitions — of  the  total 
energy  has  been  demonstrated  close  to  the  exceptional  points 
in  Ref.  [18].  In  a  silicon  platform  claims  have  been  made  that 
nonreciprocal  light  propagation  in  a  silicon  photonic  circuit  has 
been  recorded  [19].  VT  -synthetic  materials  can  exhibit  several 
intriguing  features.  These  include,  among  others,  power  oscil¬ 
lations  and  nonreciprocity  of  light  propagation  [11,16,20,21], 
nonreciprocal  Bloch  oscillations  [22],  and  unidirectional  in¬ 
visibility  [23].  More  specifically,  a  recent  paper  has  proposed 
photonic  honeycomb  lattices  with  VT  symmetry  [10].  Inter¬ 
estingly,  that  work  has  shown  that  introducing  alternating  gain 
and  loss  to  a  honeycomb  system  prohibits  VT  symmetry,  but 
adding  appropriate  strain  (direction  and  strength)  restores  the 
symmetry,  giving  rise  to  VT  -symmetric  photonic  lattices  [10]. 
Moreover,  the  mentioned  work  has  found  that  in  such  systems, 
much  higher  group  velocities  can  be  achieved  (compared 
with  non-PT-symmetry  breaking  systems),  corresponding  to 
a  tachyonic  dispersion  relation  [10].  In  the  nonlinear  domain, 
such  pseudo-Hermitian  nonreciprocal  effects  can  be  used  to 
realize  a  new  generation  of  on-chip  isolators  and  circulators 
[24] .  Other  results  within  the  framework  of  VT  optics  include 
the  realization  of  coherent  perfect  laser  absorber  [25],  spatial 
optical  switches  [26],  and  nonlinear  switching  structures  [27]. 
Despite  the  wealth  of  results  on  transport  properties  of  VT  - 
symmetric  one-dimensional  optical  structures,  the  properties 
of  high-dimensional  VT  optical  lattices  (with  the  exception 
of  few  recent  studies  [11,28])  has  remained  so  far  essentially 
unexplored. 
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Recently,  it  was  pointed  out  [29]  that  PT-symmetric 
Hamiltonians  are  a  special  case  of  pseudo-Hermitian  Hamil¬ 
tonians,  i.e.,  Hamiltonians  that  have  an  antilinear  symmetry 
[30-32].  Such  Hamiltonians  commute  with  an  antilinear 
operator  ST ,  where  S  is  a  generic  linear  operator.  The 
corresponding  Hamiltonian  TL  is  termed  generalized  VT 
symmetric.  In  a  similar  manner  as  in  the  case  of  VT  symmetry, 
one  finds  that  if  the  eigenstates  of  TL  are  also  eigenstates  of  the 
ST  operator  then  all  the  eigenvalues  of  TL  are  strictly  real  and 
the  ST  symmetry  is  said  to  be  exact.  Otherwise  the  symmetry 
is  said  to  be  broken.  An  example  case  of  a  generalized 
PT-symmetric  optical  structure  is  the  one  reported  in  [33] 
(see  also  [20]).  The  unifying  feature  of  these  systems  is  that 
they  are  built  of  a  particular  kind  of  “building  block”  which 
below  is  referred  to  as  a  dimer.  Each  dimer  in  itself  does  have 
VT  symmetry  and  it  can  be  represented  as  a  pair  of  sites 
with  assigned  energies  en ,  e  * .  The  system  is  composed  of  such 
dimers,  coupled  in  some  way.  For  an  arbitrary  choice  of  the 
site  energies  €n  and  coupling  between  the  dimers,  the  system 
as  a  whole  does  not  possess  VT  symmetry  (indeed,  such 
global  VT  symmetry  would  require  precise  relation  between 
various  en  and  coupling  symmetry  between  the  dimers).  On 
the  other  hand,  since  each  dimer  is  VT  symmetric,  with 
respect  to  its  own  center,  there  is  some  kind  of  “local”  VT 
symmetry  (which  we  shall  define  as  VdT  symmetry).  The 
main  message  of  these  papers  is  that,  VdT  symmetry  ensures  a 
robust  region  of  parameters  in  which  the  system  has  an  entirely 
real  spectrum.  Below,  we  will  use  the  term  VT  symmetry 
in  a  loose  manner  and  we  will  include  also  systems  with 
generalized  VT  symmetry. 

In  this  paper  we  investigate  beam  propagation  in  non- 
Hermitian  two-dimensional  photonic  honeycomb  lattices  with 
VT  symmetry  and  probe  for  the  possibility  of  abnormal 
diffraction.  We  find  a  type  of  conical  diffraction  that  is 
associated  with  the  spontaneously  PT-symmetry -breaking 
phase  transition  point.  Despite  the  fact  that  at  the  EP  the  Hilbert 
space  collapses,  the  emerging  cone  is  brighter  and  propagates 
with  a  transverse  velocity  that  is  controlled  by  the  gain  and 
loss  parameter  y . 

The  organization  of  the  paper  is  as  follows:  In  Sec.  II,  we 
present  the  mathematical  model.  In  Sec.  Ill  we  analyze  the 
stationary  properties  of  the  system  and  introduce  a  criterion 
based  on  the  degree  of  nonorthogonality  of  the  eigenmodes 
in  order  to  identify  the  EP  for  finite  system  sizes.  The  beam 
dynamics  is  studied  in  Sec.  IV.  First,  we  study  numerically  the 
beam  evolution  in  Sec.  IV  A,  while  the  theoretical  analysis  is 
performed  in  Sec.  IV  B.  Our  conclusions  are  given  in  Sec.  V. 


II.  PHOTONIC  LATTICE  MODEL 

We  consider  a  two-dimensional  honeycomb  photonic  lattice 
of  coupled  optical  waveguides.  Each  waveguide  supports 
only  one  mode,  while  light  is  transferred  from  waveguide 
to  waveguide  through  optical  tunneling.  A  schematic  of  the 
setup  is  shown  in  Fig.  1.  The  lattice  consist  of  two  types 
of  waveguides:  type  A  is  made  from  lossy  material  (green), 
whereas  type  B  exhibits  an  equal  amount  of  gain  (red). 
Their  arrangement  in  space  is  such  that  they  form  coupled 
A-B  dimers  with  inter-  and  intradimer  couplings  ta  and  t, 
respectively.  Such  structure,  apart  from  a  global  VT  symmetry, 
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FIG.  1 .  (Color  online)  Honeycomb  photonic  lattice  structure  with 
intradimer  coupling  t  and  interdimer  coupling  4  =  1.  Sublattice 
(lossy  waveguide)  anjfn  is  shown  by  green  (gray  in  the  black-and-white 
printed  version)  circles  while  sublattice  (gain  waveguide)  bn,m  is 
shown  by  the  red  (dark  in  the  black-and-white  printed  version)  circles. 
Each  dimer  is  distinguished  by  index  n  and  m.  The  field  is  coupled 
evanescently  between  the  waveguides. 


respects  also  another  antilinear  symmetry  (in  Ref.  [33]  we 
coined  this  VdT  symmetry)  which  is  related  to  the  local  VT 
symmetry  of  each  individual  dimer. 

Without  loss  of  generality  we  rescale  everything  in  units  of 
the  interdimer  coupling  ta  =  1 .  In  the  tight  binding  description, 
the  diffraction  dynamics  of  the  mode  electric  field  amplitude 
^n  m  —  (an^m,bn^m)T  at  the  (n,m) th  dimer  evolves  according 
to  the  following  Schrodinger-like  equation: 


i  .  T-  T  bn— i  ffl|i  T  bn— i  m-i  T  tbn  m  —  0, 

dz 


.  dbn_ 

I  — ~  +  6  bn,m  T-  Cln+\,m+\  T"  dn+\,m—l  T"  tan  m  —  0, 
dz 


(1) 


where  €  =  €o  +  iy  is  related  to  the  refractive  index  [11]. 
Without  loss  of  generality,  we  assume  60  =  0  and  y  >  0. 

It  is  useful  to  work  in  the  momentum  space.  To  this  end,  we 
write  the  field  amplitudes  ^ n,m  in  their  Fourier  representation, 
i.e., 


bn,m(T)  — 


1  f 

“  W  J_n 

1  f 

“  (2 rfj-n 


dkxdkydkx,ky(z)  exp(i[nkx  +mky]), 


(2) 


dkxdkybkx,ky(z)  enp(i[nkx  +  mky]). 


Substitution  of  Eq.  (2)  into  Eq.  (1)  leads  to  the  following  set 
of  coupled  differential  equations  in  the  momentum  space: 


dz  VMz)/  VMz)/ 

where 

Hk  =  (~ir  °<k) 

VD(k)*  iy 

and 

D(k)  =  -(t  +  2e~ikx  cos  ky);  k  =  (kx,ky). 


(3) 


(4) 

(5) 
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In  other  words,  because  of  the  translational  invariance  of  the 
system,  the  equations  of  motion  in  the  Fourier  representation 
break  up  into  2x2  blocks,  one  for  each  value  of  momentum  k. 
The  two-component  wave  functions  for  different  k  values  are 
then  decoupled,  thus  allowing  a  simple  theoretical  description 
of  the  system. 


III.  EIGENMODES  ANALYSIS 


We  start  our  analysis  with  the  study  of  the  stationary 
solutions  corresponding  to  Eq.  (3).  Substituting  the  stationary 
form 

(an,m,b„,m)T  =  exp  {-i£z){A,B)T  (6) 

in  Eq.  (3),  we  get 


The  spectrum  is  obtained  by  requesting  a  nontrivial  solution, 
i.e.,  (A,B)  ^  0.  The  corresponding  dispersion  relation  [10] 


has  the  form 

£±  =  ±V|£>(k)|2-]/2. 

(8) 

For  y  =  0  the  dispersion  relation  is 

£  =  ±\D(k)\, 

(9) 

and  we  have  two  bands  of  width  t  +  2.  There  are  three  pairs 
of  diabolic  points  (DPs): 


kj ,T  =  (  ±7 r,  zb  arccos  ^  j  , 


k  f  = 


0,  =b  (  7r  —  arccos  - 


(10) 


Expansion  of  D( k)  up  to  the  first  order  around  the  DPs  leads 
to 


£>(k)  ss  i t r]x  =f  \  —  V 


Ily¬ 


in) 


Substituting  the  above  expression  in  the  energy  dispersion 
given  by  Eq.  (8),  we  get  the  linear  relation 

£k*±[t2r12x  +  (4-t2)r12y]1/2,  (12) 

where  r]x,y  =  kx,y  -  (k0,i)X)y. 

The  standard  passive  (y  =  0)  honeycomb  lattice  (zero 
strain)  corresponds  to  t  =  1 .  In  this  case,  there  are  three  pairs 
of  DPs  at  [see  Fig.  2(a)] 


Fori  <  t  <  2  the  two  pairs  of  DPs  at  kj ,T  start  moving  toward 
each  other  while  the  pair  at  k^  moves  away  from  one  another. 
At  t  =  2  a  degeneracy  occurs,  i.e.,  kj ,T  =  (zb7r,0).  At  the 
same  time  the  dispersion  relation  £  around  kj ,T  and  kf  is 
linear  only  in  the  kx  direction  (and  quadratic  in  the  ky).  For 
t  >  2  the  two  energy  surfaces  move  away  from  each  other  and 
a  gap  between  them  is  created  [see  Fig.  2(b)].  Therefore  for 
t  ^  2  the  DPs  disappear  for  y  =  0,  and  the  conical  diffraction 
is  destroyed  [6,34]. 


FIG.  2.  (Color  online)  Energy  surfaces  for  a  honeycomb  lattice 
with  ta  =  1  and  (a)  (t,y)  =  (1,0)  where  a  DP  is  present,  (b)  (t,y)  = 
(2.5,0)  where  a  band  gap  has  been  created,  and  (c)  (t,y)  =  (2.5, 0.5) 
where  a  gain  and  loss  parameter  led  to  the  creation  of  a  EP. 

By  introducing  gain  and  loss  to  the  system  described  by 
Eq.  (1),  the  resulting  effective  Hamiltonian  which  describes 
the  paraxial  evolution  becomes  non-Hermitian.  In  fact,  for 
1  ^  t  ^  2,  any  value  of  y  results  in  complex  eigenvalues, 
i.e.,  the  system  is  in  the  broken  PT-symmetry  phase.  The 
resulting  dispersion  relation  resembles  the  dispersion  relation 
of  relativistic  particles  with  imaginary  mass,  as  was  recently 
discussed  in  Ref.  [10]. 

In  the  case  of  t  >  2  the  size  of  the  gap  between  the  two 
bands  can  be  controlled  by  manipulating  the  gain  and  loss 
parameter  y .  In  this  case,  there  is  a  y  domain,  corresponding 
to  the  exact  phase,  for  which  the  energies  are  real.  It  turns  out 
from  Eq.  (8)  that  the  line 


Yvr  =  t  —  2  (14) 

defines  the  phase  transition  from  exact  to  broken  VT  sym¬ 
metry  [10].  The  mechanism  for  this  symmetry  breaking  is  the 
crossing  between  levels,  associated  with  the  exceptional  points 

k^p  =  (0,  ±7i) 

or  (15) 

kfp  =  {±71 ,0) 

[see  Fig.  2(c)]  and  belonging  to  different  bands  [33]:  it  follows 
from  Eq.  (8),  that  when  y  =  y-pr,  the  gap  disappears  and  the 
two  (real)  levels  at  the  “inner”  band  edges  become  degenerate. 
Evaluation  of  D(kx,ky )  to  second  order  in  (r]x,r]y)  around  the 
degeneracy  points  leads  to 

D{r]x,r]y)  «b  ~(yVT  +  2i r]x  +  rj2),  rj2  =  r]2x  +  r]2,  (16) 

thus  resulting  in  the  dispersion  relation 

£  =  ±^2 y^rn 2  +  (2  ypr  +  4)^2.  (17) 

For  large  y<pr  values  (i.e.,  y<pr  2)  one  can  approximate  the 
above  equation  to  get 


£  «  ±y/2y^rj.  (18) 

This  comment  will  become  important  in  the  analysis  of  optical 
beam  propagation  discussed  in  Sec.  IV. 

Next,  we  turn  to  the  analysis  and  characterization  of  the 
biorthogonal  set  of  eigenvectors  of  our  non-Hermitian  system. 
The  target  here  is  to  identify  the  proximity  to  the  exceptional 
point  in  the  case  of  finite  Hilbert  spaces,  where  finite- size 
effects  might  play  an  important  role  in  the  analysis  of  the 
dynamics.  The  latter  do  not  respect  the  standard  (Euclidian) 
orthonormalization  condition.  Let  ( Ln\  and  \Rn)  denote  the 
left  and  right  eigenvectors  of  the  non-Hermitian  Hamiltonian 
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H  corresponding  to  the  eigenvalue  £n,  i.e., 
(Ln\H  =  (Ln\£n, 
H\Rn)  =£n\Rn)- 
The  vectors  can  be  normalized  to  satisfy 

(Ln  |  Rm )  =  Snm , 


J2\Rn)(Ln\  =  l. 

n 

Above,  J\f  is  the  dimension  of  the  Hilbert  space. 

An  observable  that  measures  the  nonorthogonality  of  the 
modes  and  can  be  used  to  identify  the  proximity  to  the  EP  in 
the  presence  of  finite-size  effects  is  the  so-called  Petermann 
factor,  which  is  defined  as  [35] 

Knm  =  (Ln\Lm)(Rm\Rn).  (22) 

We  have  studied  the  mean  (diagonal)  Petermann  factor 

-  i  £ 

k  =  j^J2k™’  (23) 

n— 1 

which  takes  the  value  of  unity  if  the  eigenfunctions  of  the 
system  are  orthogonal,  while  it  is  larger  than  unity  in  the 
opposite  case.  At  the  EP,  a  pair  of  eigenvectors  associated  with 
the  corresponding  degenerate  eigenvalues  coalesce,  leading  to 
a  collapse  of  the  Hilbert  space.  At  this  point,  the  Petermann 
factor  diverges  as  [20,36] 

K~\/\y-ypt\.  (24) 

.  This  indicates  strong  correlations  between  the  spectrum  and 
the  eigenvectors,  which  can  affect  drastically  the  dynamics,  as 
we  will  see  later. 

In  Fig.  3  we  present  our  calculations  for  K  for  different 
(t,y)  values  and  various  system  sizes.  We  find  that  as  A f 
increases,  the  divergence  approaches  the  line  y-pr  =  t  —  2, 
which  was  derived  previously  [see  Eq.  (14)]  for  the  case  of 
infinite  graphene. 


(19) 

(20) 

(21) 


IV.  DYNAMICS 

Armed  with  the  previous  knowledge  about  the  eigenmode 
properties  of  the  VT -symmetric  graphene,  we  are  now  ready 


FIG.  3.  (Color  online)  Inverse  Petermann  factor  for  various 
system  sizes  (a)  J\f  =  24,  (b)  J\f  =  168,  and  (c)  J\f  =  440.  The 
asymptotic  line  yVr  =  t  —  2  (black  solid  line)  is  approached  by  our 
numerical  data  (zero  of  K  )  for  increasing  J\f .  The  green  dotted  line 
is  plotted  to  guide  the  eye.  The  insets  show  representative  K  ’s  vs 
y  for  a  fixed  coupling  t  =  2.5. 
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to  study  beam  propagation  in  PT-symmetric  honeycomb 
lattices  at  the  vicinity  of  the  EPs.  The  question  at  hand  is 
whether  the  collapse  of  the  Hilbert  space  at  the  EP  will  affect 
the  conical  diffraction  (CD)  pattern,  and  if  yes,  what  is  the 
emerging  dynamical  picture. 

A.  Numerical  analysis 

We  first  study  wave  propagation  in  the  honeycomb  lattice 
numerically  (Fig.  4),  by  launching  a  beam  with  the  structure  of 
a  Bloch  mode  associated  with  the  EP,  multiplied  by  a  Gaussian 
envelope.  The  Bloch  modes  at  the  tip  can  be  constructed 
from  pairs  of  plane  waves  with  k  vectors  of  opposite  pairs 
of  exceptional  points.  Thus,  interfering  two  plane  waves  at 
angles  associated  with  opposite  EP  yields  the  phase  structure 
of  the  modes  from  these  points.  Multiplying  these  waves 
by  an  envelope  yields  a  superposition  of  Bloch  modes  in  a 
region  around  these  points.  Figure  4  shows  an  example  of  the 
propagation  of  a  beam  constructed  to  excite  a  Gaussian  super¬ 
position  of  Bloch  modes  around  an  EP.  The  input  beam  has  a 
bell-shape  structure,  which,  after  some  distance,  transforms 
into  the  ring-like  characteristic  of  conical  diffraction  [34]. 
From  there  on,  the  ring  is  propagating  in  the  lattice  by  keeping 
its  width  constant  while  its  radius  is  increasing  linearly  with 
distance.  The  invariance  of  the  ring  thickness  and  structure 
manifests  a  (quasi-)  linear  dispersion  relation  above  and  below 
the  EP  (see  Fig.  2);  hence,  the  diffraction  coefficient  for  wave 
packets  constructed  from  Bloch  modes  in  that  region  is  zero 
(infinite  effective  mass).  This  is  especially  interesting  because 
the  ring  itself  is  a  manifestation  of  the  dispersion  properties 
at  the  EP  itself,  where  the  diffraction  coefficient  is  infinite 
(zero-effective  mass).  As  a  result,  the  ring  forms  a  light  cone 
in  the  lattice.  The  appearance  of  CD  in  the  case  of  VT  lattices 
where  the  eigenvectors  are  nonorthogonal  and  coalesce  at  the 
EP  singularity  provides  a  clear  indication  that  the  phenomenon 
is  insensitive  to  the  eigenmodes  structure  and  that  it  depends 
only  on  the  properties  of  the  dispersion  relation. 

The  VT  -symmetric  conical  diffraction  shows  some  unique 
characteristics  with  respect  to  the  CD  found  in  the  case  of 
beam  propagation  around  DPs  for  passive  honeycomb  lattices 


FIG.  4.  (Color  online)  Propagation  of  a  Gaussian  superposition 
of  Bloch  modes  associated  with  the  vicinity  of  a  EP  at  k  =  (0,  ±tt), 
in  a  VT  -symmetric  honeycomb  lattice.  Shown  is  the  beam  intensity 
at  normalized  propagation  distances  of  (a),(d)  z  =  0;  (b),(e)z  =  10; 
(c),(f)  z  =  15,  for  yVr  =  1  and  t  =  3  (upper  panels)  and  yVr  =  2  and 
t  —  4  (lower  panels).  The  input  bell-shaped  beam  transforms  into  a 
ring-like  structure  of  light  of  a  nonvarying  thickness. 
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the  field  (an,m,bn,m)  is 

(-1  )li[z  ~  yvT(p(n,m,z,g)]  +  g 


1=1,2 


bn,m(,Z)  ^  ^  ^ 
1= 1,2 


4>(n,m,z,g )3/2 

(-l?+1[ypT<Kn,m,z,g)+  ^=|] 
<p{n,m,z,g )3/2 


(29) 


where 


FIG.  5.  (Color  online)  Transverse  velocity  of  the  spreading  ring 
vs  the  yVr  •  Numerical  simulations  approve  theoretical  predictions  of 
the  transverse  speed  of  the  CD,  which  is  vg  a  ^/yvr- 


4>{n,m,z,g)  =  [g  +  (-l)liz]2  +  n2/(2yVT+4) 

+  m2/(2y-pr).  (30) 


(i.e.,  y  =  0).  A  profound  difference  is  associated  with  the  fact 
that  now  the  transverse  speed  of  the  cone  is  increased  [10], 
and  in  fact  it  can  be  controlled  by  the  magnitude  of  the  gain 
and  loss  parameter  at  the  symmetry-breaking  point  ypr-  This 
is  shown  in  Fig.  4  where  we  compare  the  spreading  of  a  CD 
for  two  different  y^r  values. 

In  Fig.  5,  we  report  in  a  double  logarithmic  plot,  our 
numerical  measurements  for  the  transverse  velocity  of  the 
spreading  ring  for  various  y^r  values.  The  best  linear  fitting 
to  the  numerical  data  shows  that  transverse  speed  of  the  cone 
is  proportional  to  ^/yrr-  This  behavior  can  be  understood 
qualitatively  by  realizing  that  the  group  velocity  near  the  EPs  is 

vg  =  dS/dk  ~  V2ypr  (25) 

[see  Eqs.(17)  and  (18)].  At  the  same  time  we  find  that  the 
resulting  VT  cone  is  brighter  with  respect  to  the  one  found 
in  passive  lattices  [5]  (i.e.,  the  field  intensity  of  the  conical 
wavefront  is  larger). 

B.  Theoretical  considerations 

It  is  possible  to  gain  valuable  insight  into  the  features  of 
PT-conical  diffraction  by  considering  the  field  evolution  in 
the  momentum  space.  We  consider  for  simplicity  an  initial 
distribution  [<3k(0),£k(0)]r  that  is  symmetric  around  the  EP 
while  it  decays  exponentially  away  from  it.  Specifically  we 
assume 

[ak(0),M0)]r  =  e~g^H(l,0)T.  (26) 

Next,  we  calculate  the  evolution  matrix  U  =  e_?HkZ  where  Hk 
is  given  by  Eq.  (4).  After  a  straightforward  algebra  and  using 
the  fact  that 

Hk2  =  S2  x  1,  (27) 


Although  our  simplified  calculations  are  not  able  to  capture 
all  the  features  of  the  propagating  cone,  the  above  expression 
encompasses  the  main  characteristics  of  the  conical  diffraction 
that  we  have  observed  in  our  numerical  simulations.  At  z  =  0, 
Eq.  (29)  resembles  a  Lorentzian,  which  slowly  transforms  into 
a  ring  of  light,  whose  radius  expands  linearly  with  z  with 
velocity  *f2ypr,  while  its  thickness  remains  unchanged.  At 
the  same  time  the  field  intensity  on  the  ring  in  the  case  of 
VT  -symmetric  lattices  is  brighter  than  the  one  corresponding 
to  passive  honeycomb  lattices  (i.e.,  l/z2  vs  1/z4  behavior, 
respectively). 


V.  SUMMARY  AND  CONCLUDING  REMARKS 

We  studied  numerically  and  analytically  the  propagation 
of  waves  in  PT-honeycomb  photonic  lattices,  demonstrating 
the  existence  of  conical  diffraction  arising  solely  from  the 
presence  of  a  spontaneous  PT-symmetry-breaking  phase 
transition  point.  In  spite  of  the  fact  that  the  eigenvectors  are 
nonorthogonal  and  there  is  a  collapse  of  the  Hilbert  space  at 
the  EP,  the  emerging  cone  is  brighter  and  moves  faster  than 
the  corresponding  one  of  the  passive  structure.  Although,  the 
realization  of  such  photonic  structures  is  currently  a  challeng¬ 
ing  task,  active  electronic  circuits,  like  the  one  proposed  in 
Ref.  [18],  can  be  proven  useful  alternatives  that  might  allow 
us  to  investigate  experimentally  wave  propagation  in  extended 
PT-symmetric  lattices. 

These  findings  raise  several  intriguing  questions.  For 
example,  how  does  nonlinearity  affect  VT  -symmetric  conical 
diffraction?  What  is  the  effect  of  disorder  [37]?  Is  this  behavior 
generic  for  any  system  at  the  spontaneously  VT  -symmetry- 
breaking  point?  These  intriguing  questions  are  universal 
and  relate  to  any  field  in  which  waves  can  propagate  in  a 
periodic  potential.  It  is  expected  that  in  active  metamaterials 
such  phenomena  will  be  present  and  they  may  have  specific 
technological  importance. 


where  1  is  the  unity  matrix,  we  get 

U  =  cos(z|£|)l  -  i[sin(z\S \)/\£ |]Hk.  (28) 

Equation  (28)  is  the  starting  point  of  our  analysis.  Substituting 
Eqs.  (4),  (16),  and  (17),  we  find  that  the  evolving  amplitude  of 
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